% Script author: Ran Wu 

%STEP 1: calculate all the variables for all the years at sector level
%STEP 2: calculate variables at country level, country-sector level, world level and world-sector level
%STEP 3: repeat the analysis, this time excluding the four primary sectors

%STEP 1: calculate all the variables
%the following 11 variables are calculated:
%(1)grossimport:gross imports 
%(2)grossexport:gross exports
%(3)eefe:domestic emissions embodied in foreign final demand
%(4)eefi:foreign emissions embodied in domestic final demand
%(5)pbe:production-based emissions (without household emissions)
%(6)cbe:consumption-based emissions (without household emissions)
%(7)eege:emissions embodied in gross exports 
%(8)eegi:emissions embodied in gross imports 
%(9)eagi:emissions avoided by gross imports 
%(10)beagt:balance of avoided emissions (net onshoring)
%(11)heege:scale-adjusted emissions embodied in gross exports

%result of 2000
load('wiot00.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity

finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end

interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end

tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export

eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 

Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end

eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export

teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import

gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector

grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector

eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end

beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade

scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end

excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)
aresult00=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult00','aresult00');%result of 2000

%the following steps will calculate results of different years one by one
%result of 2001
load('wiot01.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult01=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult01','aresult01'); %result of 2001

%result of 2002
load('wiot02.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult02=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult02','aresult02'); %result of 2002

%result of 2003
load('wiot03.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult03=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult03','aresult03'); %result of 2003

%result of 2004
load('wiot04.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult04=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult04','aresult04'); %result of 2004

%result of 2005
load('wiot05.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult05=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult05','aresult05'); %result of 2005

%result of 2006
load('wiot06.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult06=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult06','aresult06'); %result of 2006

%result of 2007
load('wiot07.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult07=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult07','aresult07'); %result of 2007

%result of 2008
load('wiot08.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult08=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult08','aresult08'); %result of 2008

%result of 2009
load('wiot09.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult09=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult09','aresult09'); %result of 2009

%result of 2010
load('wiot10.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult10=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult10','aresult10'); %result of 2010

%result of 2011
load('wiot11.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult11=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult11','aresult11'); %result of 2011

%result of 2012
load('wiot12.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult12=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult12','aresult12'); %result of 2012

%result of 2013
load('wiot13.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult13=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult13','aresult13'); %result of 2013

%result of 2014
load('wiot14.mat')
m=56;
n=44;
Z=wiot(1:m*n,1:m*n);%direct consumption coefficient matrix Z
X=wiot(1:m*n,m*n+n*5+1);%total output vector X
Y=wiot(1:m*n,m*n+1:m*n+n*5);%final demand matrix Y m*n*n*5
for i=1:m*n
    A(:,i)=Z(:,i)/X(i);
end
A(isnan(A))=0;
I=eye(m*n);
B=inv(I-A);%Leontief inverse matrix B
%during pretreatment process, emissions data(in million tonnes) have been merged with WIOT input-output tables, as shown in column2686
emissionmt=wiot(1:m*n,m*n+n*5+2);%emission vector at country-sector level(in mt)
pbe=emissionmt;%production-based emissions (without household emissions)
intensity=emissionmt./X;%emission intensity of total output
intensity(isnan(intensity))=0;
intensity(intensity==inf)=0;
V=diag(intensity);%diagonalization of emission intensity
 
finaldemand=Y;%rename matrix Y with finaldemand
finalsim=zeros(m*n,n);%transfer final demand matrix from (m*n)*(n*5) to (m*n)*n 
for i=0:n-1
    finalsim(:,1+i)=sum(finaldemand(:,(1+5*i):(5+5*i)),2);
end
 
interinput=Z;%rename intermediate input matrix Z with interinput
interinputsim=zeros(m*n,n);%transfer intermediate input matrix from m*nxm*n to m*n*n
for i=0:n-1
    interinputsim(:,1+i)=sum(interinput(:,(1+m*i):(m+m*i)),2);
end
 
tradeflow=finalsim+interinputsim;%trade flow matrix
eefinal=V*B*finalsim;%matrix of emission embodied in final trade flow
for i=0:n-1
    eefinal((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eefe=sum(eefinal,2);%column vector-EEFE:emission embodied in final export
 
eefimatrix=zeros(m,n);%matrix of emission embodied in final import 
for i=0:n-1
    matrix1=eefinal((1+m*i):(m+m*i),:);
    eefimatrix=eefimatrix+matrix1;
    i=i+1;
end
eefi=eefimatrix(:);%column vector-EEFI:emission embodied in final import 
 
Bdomestic=zeros(m*n,m*n);%keep domestic part only in Loentif inverse matrix B
for i=0:n-1
    Bdomestic((1+m*i):(m+m*i),(1+m*i):(m+m*i))=B((1+m*i):(m+m*i),(1+m*i):(m+m*i));
end
 
eegross=V*Bdomestic*tradeflow;%matrix of domestic emissions embodied in trade flow
for i=0:n-1
    eegross((1+m*i):(m+m*i),i+1)=0;%set domestic part as 0   
end
eege=sum(eegross,2);%column vector-EEGE:domestic emissions embodied in gross export
 
teefimatrix=zeros(m,n);%matrix of foreign emissions embodied in gross import
for i=0:n-1
    matrix3=eegross((1+m*i):(m+m*i),:);
    teefimatrix=teefimatrix+matrix3;
    i=i+1;
end
eegi=teefimatrix(:);%column vector-EEGI:foreign emissions embodied in gross import
 
gross=tradeflow;%generate matrix gross=tradeflow
for i=0:n-1
    gross((1+m*i):(m+m*i),1+i)=0; %set domestic part as 0      
end
grossexport=sum(gross,2);%gross export vector
 
grossimportmatrix=zeros(m,n);
for i=0:n-1
    gross1=gross((1+m*i):(m+m*i),:);
    grossimportmatrix=grossimportmatrix+gross1;
    i=i+1;
end
grossimport=grossimportmatrix(:);%gross import vector
 
eagi=V*Bdomestic*grossimport;%domestic emissions avoid by gross import
for i=1:m*n
    if grossimport(i)~=0 && intensity(i)==0
        eagi(i)=eegi(i);
    else
        eagi(i)=eagi(i);
    end
end
 
beagt=eege-eagi;%BEAGT-balance of emissions avoided by trade
 
scaleim=zeros(m*n,1);%hypothetical EEGE using the actual export composition but the scale of imports
for i=0:n-1
    scaleim((1+m*i):(m+m*i),1)=sum(grossimport((1+m*i):(m+m*i),1),1);
end
scaleex=zeros(m*n,1); 
for i=0:n-1
    scaleex((1+m*i):(m+m*i),1)=sum(grossexport((1+m*i):(m+m*i),1),1);
end
 
excomp=grossexport./scaleex;%actual export composition
hgrossexport=excomp.*scaleim;%hypothetical gross export
heege=V*Bdomestic*hgrossexport;%scale-adjusted EEGE
cbe=pbe-eefe+eefi;%consumption-based emissions (without household emissions)

aresult14=[grossimport grossexport eefe eefi pbe cbe eege eegi eagi beagt heege];
save('aresult14','aresult14'); %result of 2014
%results of all the years have been calculated till this step

%STEP 2: calculate variables at country level, country-sector level, world level and world-sector level
%combine the results from 2000 to 2014
clear
load('aresult00.mat')
load('aresult01.mat')
load('aresult02.mat')
load('aresult03.mat')
load('aresult04.mat')
load('aresult05.mat')
load('aresult06.mat')
load('aresult07.mat')
load('aresult08.mat')
load('aresult09.mat')
load('aresult10.mat')
load('aresult11.mat')
load('aresult12.mat')
load('aresult13.mat')
load('aresult14.mat')
aresult=[aresult00;aresult01;aresult02;aresult03;aresult04;aresult05;aresult06;aresult07;aresult08;aresult09;aresult10;aresult11;aresult12;aresult13;aresult14];
save('aresult','aresult')
columns1={'grossimport','grossexport', 'eefe', 'eefi','pbe', 'cbe','eege','eegi', 'eagi', 'beagt', 'heege'};
aresultcell=num2cell(aresult);
columnsaresult=[columns1;aresultcell];
xlswrite('io_result_all.xlsx',columnsaresult,'sheet2');
%the following code is to merge results of step 1
%results at country level
load('householdmt.mat')
m=56;
n=44;
acountry=zeros(n*15,11);%15years,11variables
for i=0:n*15-1
    acountry(i+1,:)=sum(aresult((1+m*i):(m+m*i),:));
end
acountry=[acountry householdmt];
save('acountry','acountry')
columns2={'grossimport','grossexport', 'eefe', 'eefi','pbe', 'cbe','eege','eegi', 'eagi', 'beagt', 'heege','household'};
acountrycell=num2cell(acountry);
columnsacountry=[columns2;acountrycell];
xlswrite('io_result_all.xlsx',columnsacountry,'sheet3');
%results at world level
aworld=zeros(15,12);
for i=0:14 %14=15-1
    aworld(i+1,:)=sum(acountry((1+n*i):(n+n*i),:));
end
save('aworld','aworld')
aworldcell=num2cell(aworld);
columnsaworld=[columns2;aworldcell];
xlswrite('io_result_all.xlsx',columnsaworld,'sheet4');

%separate results of all the years to calculate world sector level result
part00=aresult(1:m*n,:);
part01=aresult(m*n+1:m*n*2,:);
part02=aresult(m*n*2+1:m*n*3,:);
part03=aresult(m*n*3+1:m*n*4,:);
part04=aresult(m*n*4+1:m*n*5,:);
part05=aresult(m*n*5+1:m*n*6,:);
part06=aresult(m*n*6+1:m*n*7,:);
part07=aresult(m*n*7+1:m*n*8,:);
part08=aresult(m*n*8+1:m*n*9,:);
part09=aresult(m*n*9+1:m*n*10,:);
part10=aresult(m*n*10+1:m*n*11,:);
part11=aresult(m*n*11+1:m*n*12,:);
part12=aresult(m*n*12+1:m*n*13,:);
part13=aresult(m*n*13+1:m*n*14,:);
part14=aresult(m*n*14+1:m*n*15,:);

%world-sector level results of 2000 to 2014
world00=zeros(m,11);
for i=0:n-1
    p00=part00((1+m*i):(m+m*i),:);
    world00=world00+p00;
    i=i+1;
end

world01=zeros(m,11);
for i=0:n-1
    p01=part01((1+m*i):(m+m*i),:);
    world01=world01+p01;
    i=i+1;
end

world02=zeros(m,11);
for i=0:n-1
    p02=part02((1+m*i):(m+m*i),:);
    world02=world02+p02;
    i=i+1;
end

world03=zeros(m,11);
for i=0:n-1
    p03=part03((1+m*i):(m+m*i),:);
    world03=world03+p03;
    i=i+1;
end
 
world04=zeros(m,11);
for i=0:n-1
    p04=part04((1+m*i):(m+m*i),:);
    world04=world04+p04;
    i=i+1;
end

world05=zeros(m,11);
for i=0:n-1
    p05=part05((1+m*i):(m+m*i),:);
    world05=world05+p05;
    i=i+1;
end

world06=zeros(m,11);
for i=0:n-1
    p06=part06((1+m*i):(m+m*i),:);
    world06=world06+p06;
    i=i+1;
end

world07=zeros(m,11);
for i=0:n-1
    p07=part07((1+m*i):(m+m*i),:);
    world07=world07+p07;
    i=i+1;
end
 
world08=zeros(m,11);
for i=0:n-1
    p08=part08((1+m*i):(m+m*i),:);
    world08=world08+p08;
    i=i+1;
end
 
world09=zeros(m,11);
for i=0:n-1
    p09=part09((1+m*i):(m+m*i),:);
    world09=world09+p09;
    i=i+1;
end

world10=zeros(m,11);
for i=0:n-1
    p10=part10((1+m*i):(m+m*i),:);
    world10=world10+p10;
    i=i+1;
end

world11=zeros(m,11);
for i=0:n-1
    p11=part11((1+m*i):(m+m*i),:);
    world11=world11+p11;
    i=i+1;
end

world12=zeros(m,11);
for i=0:n-1
    p12=part12((1+m*i):(m+m*i),:);
    world12=world12+p12;
    i=i+1;
end

world13=zeros(m,11);
for i=0:n-1
    p13=part13((1+m*i):(m+m*i),:);
    world13=world13+p13;
    i=i+1;
end

world14=zeros(m,11);
for i=0:n-1
    p14=part14((1+m*i):(m+m*i),:);
    world14=world14+p14;
    i=i+1;
end

%combine the world-sector level results from 2000 to 2014
aworldsector=[world00;world01;world02;world03;world04;world05;world06;world07;world08;world09;world10;world11;world12;world13;world14];
save('aworldsector','aworldsector')
aworldsectorcell=num2cell(aworldsector);
columnsaworldsector=[columns1;aworldsectorcell];
xlswrite('io_result_all.xlsx',columnsaworldsector,'sheet5');

%STEP 3:exclude four primary industries and get results at country level and world level
clear
load('aresult.mat')
load('householdmt.mat')
m=56;
n=44;
zero=zeros(4,11);%4 primary industries, 11 variables
aresult1=aresult;
for i=0:n*15-1
   aresult1(1+m*i:4+m*i,:)= zero(1:4,:);
end
save('aresult1','aresult1')
columns1={'grossimport','grossexport', 'eefe', 'eefi','pbe', 'cbe','eege','eegi', 'eagi', 'beagt', 'heege'};
aresultcell1=num2cell(aresult1);
columnsaresult1=[columns1;aresultcell1];
xlswrite('io_result_wo_primary.xlsx',columnsaresult1,'sheet2');

acountry1=zeros(n*15,11);
for i=0:n*15-1
    acountry1(i+1,:)=sum(aresult1((1+m*i):(m+m*i),:));
end
acountry1=[acountry1 householdmt];
save('acountry1','acountry1')
columns2={'grossimport','grossexport', 'eefe', 'eefi','pbe', 'cbe','eege','eegi', 'eagi', 'beagt', 'heege','household'};
acountrycell1=num2cell(acountry1);
columnsacountry1=[columns2;acountrycell1];
xlswrite('io_result_wo_primary.xlsx',columnsacountry1,'sheet3');

aworld1=zeros(15,12);
for i=0:14
    aworld1(i+1,:)=sum(acountry1((1+n*i):(n+n*i),:));
end
save('aworld1','aworld1')
aworldcell1=num2cell(aworld1);
columnsaworld1=[columns2;aworldcell1];
xlswrite('io_result_wo_primary.xlsx',columnsaworld1,'sheet4');

%separate results of all the years to calculate world sector level result
part00=aresult1(1:m*n,:);
part01=aresult1(m*n+1:m*n*2,:);
part02=aresult1(m*n*2+1:m*n*3,:);
part03=aresult1(m*n*3+1:m*n*4,:);
part04=aresult1(m*n*4+1:m*n*5,:);
part05=aresult1(m*n*5+1:m*n*6,:);
part06=aresult1(m*n*6+1:m*n*7,:);
part07=aresult1(m*n*7+1:m*n*8,:);
part08=aresult1(m*n*8+1:m*n*9,:);
part09=aresult1(m*n*9+1:m*n*10,:);
part10=aresult1(m*n*10+1:m*n*11,:);
part11=aresult1(m*n*11+1:m*n*12,:);
part12=aresult1(m*n*12+1:m*n*13,:);
part13=aresult1(m*n*13+1:m*n*14,:);
part14=aresult1(m*n*14+1:m*n*15,:);

%world-sector level results of 2000 to 2014
world00=zeros(m,11);
for i=0:n-1
    p00=part00((1+m*i):(m+m*i),:);
    world00=world00+p00;
    i=i+1;
end

world01=zeros(m,11);
for i=0:n-1
    p01=part01((1+m*i):(m+m*i),:);
    world01=world01+p01;
    i=i+1;
end

world02=zeros(m,11);
for i=0:n-1
    p02=part02((1+m*i):(m+m*i),:);
    world02=world02+p02;
    i=i+1;
end

world03=zeros(m,11);
for i=0:n-1
    p03=part03((1+m*i):(m+m*i),:);
    world03=world03+p03;
    i=i+1;
end
 
world04=zeros(m,11);
for i=0:n-1
    p04=part04((1+m*i):(m+m*i),:);
    world04=world04+p04;
    i=i+1;
end

world05=zeros(m,11);
for i=0:n-1
    p05=part05((1+m*i):(m+m*i),:);
    world05=world05+p05;
    i=i+1;
end

world06=zeros(m,11);
for i=0:n-1
    p06=part06((1+m*i):(m+m*i),:);
    world06=world06+p06;
    i=i+1;
end

world07=zeros(m,11);
for i=0:n-1
    p07=part07((1+m*i):(m+m*i),:);
    world07=world07+p07;
    i=i+1;
end
 
world08=zeros(m,11);
for i=0:n-1
    p08=part08((1+m*i):(m+m*i),:);
    world08=world08+p08;
    i=i+1;
end
 
world09=zeros(m,11);
for i=0:n-1
    p09=part09((1+m*i):(m+m*i),:);
    world09=world09+p09;
    i=i+1;
end

world10=zeros(m,11);
for i=0:n-1
    p10=part10((1+m*i):(m+m*i),:);
    world10=world10+p10;
    i=i+1;
end

world11=zeros(m,11);
for i=0:n-1
    p11=part11((1+m*i):(m+m*i),:);
    world11=world11+p11;
    i=i+1;
end

world12=zeros(m,11);
for i=0:n-1
    p12=part12((1+m*i):(m+m*i),:);
    world12=world12+p12;
    i=i+1;
end

world13=zeros(m,11);
for i=0:n-1
    p13=part13((1+m*i):(m+m*i),:);
    world13=world13+p13;
    i=i+1;
end

world14=zeros(m,11);
for i=0:n-1
    p14=part14((1+m*i):(m+m*i),:);
    world14=world14+p14;
    i=i+1;
end
%combine world-sector level results from 2000 to 2014 (without 4 primary industries)
aworldsector1=[world00;world01;world02;world03;world04;world05;world06;world07;world08;world09;world10;world11;world12;world13;world14];
save('aworldsector1','aworldsector1')
aworldsectorcell1=num2cell(aworldsector1);
columnsaworldsector1=[columns1;aworldsectorcell1];
xlswrite('io_result_wo_primary.xlsx',columnsaworldsector1,'sheet5');

%aresult: results at country-sector level (including all the industries)
%aresult1: results at country-sector level  (without 4 primary industries)
%acountry: results at country level (including all the industries)
%acountry1: results at country level (without 4 primary industries)
%aworldsector: results at world-sector level (including all the industries)
%aworldsector1: results at world-sector level (without 4 primary industries)
%aworld: results at world level (including all the industries)
%aworld1: results at world level (without 4 primary industries) 

%Numbers in A1
clear
m=56;
n=44;
load('acountry1')
load('householdmt')
load('acountry')
pbe=acountry(:,5)+householdmt;%production-based emissions at country level
neton1=acountry1(:,7)-acountry1(:,9);%net onshoring=eege-eagi (without primary industries)
share1=neton1./pbe;
share1=reshape(share1,n,15);
neton1=reshape(neton1,n,15);
avshare1=sum(share1,2)./15;
avneton1=sum(neton1,2)./15;
neton=acountry(:,7)-acountry(:,9);%net onshoring=eege-eagi (including all the industries)
share=neton./pbe;
share=reshape(share,44,15);
neton=reshape(neton,44,15);
avshare=sum(share,2)./15;
avneton=sum(neton,2)./15;
result1=avshare1-avshare;
result2=avneton1-avneton;
appendix1=[result2 avshare1 avshare];
save('appendix1','appendix1');
columns={'mean_deviations','neton_pbe_wo', 'neton_pbe_all'};
appendixcell1=num2cell(appendix1);
columnsappendix1=[columns;appendixcell1];
xlswrite('appendix.xlsx',columnsappendix1,'sheet1');

%Results in Table 4: Trade-adjusted PBE (Appendix)
clear
m=56;
n=44;
load('acountry1')
load('acountry')
load('householdmt')
pbe=acountry(:,5)+householdmt;%production-based emissions at country level
neton=acountry1(:,7)-acountry1(:,9);%net onshoring=eege-eagi
tpbe=pbe-neton;%trade-adjusted production-based emissions
column1=tpbe(1:n,:);%results of 2000
column2=tpbe(n*14+1:n*15,:);%results of 2014
column3=column2-column1;
column4=column3./column1;
table4=[column1 column2 column3 column4];
save('table4','table4');
columns={'2000','2014', 'der_mt', 'der_percent'};
tablecell4=num2cell(table4);
columnstable4=[columns;tablecell4];
xlswrite('appendix.xlsx',columnstable4,'sheet2');

%Results in Table 8: Net onshoring (Appendix)
clear
m=56;
n=44;
load('acountry1')
load('householdmt')
load('acountry')
pbe=acountry(:,5)+householdmt;%production-based emissions at country level
neton=acountry1(:,7)-acountry1(:,9);%net onshoring=eege-eagi
share=neton./pbe;
column1=share(1:n,:);%results of 2000
column2=share(n*7+1:n*8,:);%results of 2007
column3=share(n*14+1:n*15,:);%results of 2014
column4=neton(1:n,:);%results of 2000
column5=neton(n*7+1:n*8,:);%results of 2007
column6=neton(n*14+1:n*15,:);%results of 2014
table8=[column1 column2 column3 column4 column5 column6];
save('table8','table8');
columns={'neton_pbe2000_wo','neton_pbe2007_wo','neton_pbe2014_wo', 'neton2000_wo','neton2007_wo','neton2014_wo'};
tablecell8=num2cell(table8);
columnstable8=[columns;tablecell8];
xlswrite('appendix.xlsx',columnstable8,'sheet3');

%Results in Table 9: Net onshoring including the primary industries (Appendix)
clear
m=56;
n=44;
load('acountry')
load('householdmt')
pbe=acountry(:,5)+householdmt;%production-based emissions at country level
neton=acountry(:,7)-acountry(:,9);%net onshoring=eege-eagi
share=neton./pbe;
column1=share(1:n,:);%results of 2000
column2=share(n*7+1:n*8,:);%results of 2007
column3=share(n*14+1:n*15,:);%results of 2014
column4=neton(1:n,:);%results of 2000
column5=neton(n*7+1:n*8,:);%results of 2007
column6=neton(n*14+1:n*15,:);%results of 2014
table9=[column1 column2 column3 column4 column5 column6];
save('table9','table9');
columns={'neton_pbe2000_all','neton_pbe2007_all','neton_pbe2014_all', 'neton2000_all','neton2007_all','neton2014_all'};
tablecell9=num2cell(table9);
columnstable9=[columns;tablecell9];
xlswrite('appendix.xlsx',columnstable9,'sheet4');

%Results in Table 10: Scale-adjusted net onshoring (Appendix)
clear
m=56;
n=44;
load('acountry1')
load('householdmt')
load('acountry')
pbe=acountry(:,5)+householdmt;%production-based emissions at country level
scaleneton=acountry1(:,11)-acountry1(:,9);%net onshoring=heege-eagi
share=scaleneton./pbe;
column1=share(1:n,:);%results of 2000
column2=share(n*7+1:n*8,:);%results of 2007
column3=share(n*14+1:n*15,:);%results of 2014
column4=scaleneton(1:n,:);%results of 2000
column5=scaleneton(n*7+1:n*8,:);%results of 2007
column6=scaleneton(n*14+1:n*15,:);%results of 2014
table10=[column1 column2 column3 column4 column5 column6];
save('table10','table10');
columns={'adjneton_pbe2000','adjneton_pbe2007','adjneton_pbe2014', 'adjneton2000','adjneton2007','adjneton2014'};
tablecell10=num2cell(table10);
columnstable10=[columns;tablecell10];
xlswrite('appendix.xlsx',columnstable10,'sheet5');

%Results in Table 11: BEET (Appendix)
clear
m=56;
n=44;
load('acountry1')
load('householdmt')
load('acountry')
beet=acountry1(:,5)-acountry1(:,6);
pbe=acountry(:,5)+householdmt;%production-based emissions at country level
share=beet./pbe;
column1=share(1:n,:);%results of 2000
column2=share(n*7+1:n*8,:);%results of 2007
column3=share(n*14+1:n*15,:);%results of 2014
column4=beet(1:n,:);%results of 2000
column5=beet(n*7+1:n*8,:);%results of 2007
column6=beet(n*14+1:n*15,:);%results of 2014
table11=[column1 column2 column3 column4 column5 column6];
save('table11','table11');
columns={'beet_pbe2000','beet_pbe2007','beet_pbe2014', 'beet2000','beet2007','beet2014'};
tablecell11=num2cell(table11);
columnstable11=[columns;tablecell11];
xlswrite('appendix.xlsx',columnstable11,'sheet6');




